AD-A147  719 


UNCLASSIFIED 


A  MOVING  FINITE  ELEMENT  METHOD  FOR  TIME  DEPENDENT 
PARTIAL  DIFFERENTIAL  EQ.  .  (U)  RENSSELAER  POLVTECHNIC 
INST  TROV  NY  DEPT  OF  MATHEMATICAL  SCIE.  . 

S  ADJERID  ET  AL.  SEP  84  AFOSR-TR-84-0946  F/G  12/1 


DTK  FILE  COPY 

I  Hi  a  Hi  g 


SFOSR-TR-  8  4  -0946 


A  MOVING  FINITE  ELEMENT  METHOD  FOR  TIME  DEPENDENT  PARTIAL 
DIFFERENTIAL  EQUATIONS  WITH  ERROR  ESTIMATION  AND  REFINEMENT* 


Slimane  Adjerid**  and  Joseph  E.  Flaherty*** 

Dedicated  in  memory  of  Richard  C.  DiPrima 
Abstract 


We  discuss  a  moving  finite  element  method  for  solving  vector  systems  of 
time  dependent  partial  differential  equations  in  one  space  dimension.  The  mesh 
is  moved  so  as  to  equidistribute  the  spatial  component  of  the  discretization 
error  in  H1 ,  We  present  a  method  of  estimating  this  error  by  using 
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or  the  temporal  integration  of  the  solution,  the  error  estimate,  and  the  mesh 
otion.  Computational  results  using  a  code  based  on  our  method  are  presented 
or  several  examples. 


Approved  for  public  release  f 
distribution  unlimited. 


•The  authors  were  partially  supported  by  the  U.  S.  Air  Force  Office  of 
Scientific  Research,  Air  Force  Systems  Command,  USAF,  under  Grant  Number 
AFOSR  80-0192  and  the  U.  S.  Army  Research  Office  under  Contract  Number 
DAAG29-82-K-0197.  The  first  author  also  received  support  from  the 
government  of  Algeria. 

•♦Department  of  Mathematical  Sciences,  Rensselaer  Polytechnic  Institute, 
Troy,  NY  12181.  This  work  will  be  used  in  partial  fulfillment  of  his  Ph.D. 
requirements  at  the  Rensselaer  Polytechnic  Institute. 

•••Department  of  Computer  Science,  Rensselaer  Polytechnic  Institute,  Troy, 

NY  12181 . 

84  11  16  022 


SECURITY  CLASSlF ICATIQN  OF  THIS  PAGE 


REPORT  DOCUMENTATION  PAGE 


la  REPORT  SECURITY  CLASSIPICAT  ION 

UNCLASSIFIED 

lb  RESTRICTIVE  MARKINGS 

3a.  SECuRitv  classification  authority 

3.  OISTRIBUTiON/A  VAilABilITY  OP  REPORT 

Approved  for  public  release;  distribution 

2b  DECLASSIPICATION/DOWNGRAQING  SCHEDULE 

unlimited.  p, 

4  PERFORMING  ORGANIZATION  REPORT  NUMBER(S) 

5  MONITORING  ORGANIZATION  REPORT  NUMBER. Si 

AFOSR-TR-  3  4.094  6 

ba  NAME  OP  PERFORMING  ORGANIZATION 

Rensselaer  Polytechnic 
Institute 

6o.  office  symbol 

(If  applicable  t 

7a  NAME  OF  MONITORING  ORGANIZATION 

Air  Force  Office  of  Scientific  Research 

6c  AOO«ESS  »Cifv  State  and /!P  Code  > 

Department  of  Mathematical  Sciences 

Troy  NY  12181 

* 

7b  ADDRESS  rCifv.  State  and  /lb 

Directorate  of  Mathematical  and  Information 
Sciences,  Boiling  AFE  DC  20332 

•».  NAME  Of  FUNOiNG/SPONSOAlNC 
ORGANIZATION 

AFOS? 

8b  OFFICE  SYMBOL 
(If  applicable  > 

NM 

9  PROCUREMENT  INSTRUMENT  IDENTIFICATION  NUMB E  A 

AF0SR-80-0192 

8c  AOORESS  •  Ci f> .  State  and  /IP  Code > 


line  AF3  DC  2C322 


10  SOURCE  OF  FUNO'NG  NOS 


11  t*Tlc  Include  Security  Classification  <  * 

A  MOVING  FINITE  ELEMENT  METHOD  FOR  TIME  DEPENDENT  PARTIAL  DIFFERENTIAL  EQUATIONS  V/ITH 


12  personal  AUTHORS)  /ERROR  ESTIMATION  AND  REFINEMENT _ 

Slimane  Adierid  and  Joseph  E.  Flaherty* 


14  DATE  OF  REPORT  •  Yr  .  .Vfo  .Day  IS  PAGE  CCuN* 

SEP  84  27 


13a.  T>PE  OP  REPORT  130  Time  COVERED 

Technical  prom _  to _ 


16  Supplementary  notation 

*Dept  of  Computer  Science,  RPI,  Troy  NY  12181. 


CCS A  * .  CODES  I  '8  SUBJECT  TERMS  Continue  on  ;/  nec«uar*.  ana  aenu'y  9\  biock  num'to" 

Moving  finite  element  method;  adaptive  methods;  time 
dependent  partial  differential  equations;  error  estimation; 
refinement . 


ABSTRACT  Continue  on  *rver$e  .  /  nice  Mary  and  identity  by  bloc*  numbe' 

The  authors  discuss  a  moving  finite  element  method  for  solving  vector  systems  of  time 


dependent  partial  differential  equations  in  one  space  dimension.  The  mesh  is  moved  so  as 
to  equidistribute  the  spatial  component  of  the  discretization  error  in  H  .  They  present  a 
method  of  estimating  this  error  by  using  p-hierarchic  finite  elements.  The  error  estimate 
is  also  used  in  an  adaptive  mesh  refinement  procedure  to  give  an  algorithm  that  combines 
mesh  movement  and  refinement. 


The  authors  discretize  the  partial  differential  equations  in  space  using  a  Galcrkin 
procedure  with  piecewise  linear  elements  to  approximate  the  solution  and  quadratic  element; 
to  estimate  the  error.  A  system  of  ordinary  differential  equations  for  mesh  (GONTTNUED' 


jc  distribution  a  ,  a.,. A8' lit  *  op  *as*R»c- 

vNCwASSlP  ED  uNw  M  ■(:  n  SAME  AS  a»*  _  f  C  -SEAS  ~ 


j:«  nave  z*  respgns  bn  . 

CRT  John  t.  Thomas,  J.r 


21  ABSTRAC”  SECun.T  .  . -ASSiP  ICATION 


I 


ffCUPITV  CLASSI*'CATlON  OF  THIS  FACE 


JTEM  19 ,  .ABSTRACT,  CONTINUED ^/velocities  are  used  to  control  element  motions.  The 
authors  use  existing  software  for  stiff  ordinary  differential  equations  for  the  temporal 
integration  of  the  solution,  the  error  estimate,  and  the  mesh  motion.  Computational 
results  using  a  code  based  on  this  method  are  presented  for  several  examples.,^- 


UNCLASSIFIED 


t 


1.  introduction 


AIR  FORCS  OFFICE  0?  EOTE!TTTFI?  FE?F*RCI  f AFSC) 

NOTICE  0:’  CFW*  n  r^TTC 

This  •«  ••  -  -  r<3  is 

appro-:- '  '  ;  •  ..vV  130-12. 

Di3tri’  -.-ui-itsd. 

MATTHB#  J.  XCRVin 

Chief,  teohuioal  Inf omation Division 


✓ 


Many  technological  situations  involve  the  rapid  formation,  evolution, 
propagation,  and  disintegration  of  small  scale  structures.  Some  examples  are 
shock  waves,  shear  layers  in  laminar  and  turbulent  flows,  phase  boundaries 
during  nonequilibrium  thermal  processes,  and  classical  boundary  layers.  With 
increasing  complexity  of  the  physical  problem,  there  is  an  increasing  need  for 
reliable  and  robust  software  tools  to  accurately  and  efficiently  describe  the 
phenomena.  Adaptive  techniques  automatically  change  and  evolve  with  the 
solution  and  are  thus  good  candidates  for  providing  the  computational  methods 
and  codes  necessary  to  solve  some  of  these  difficult  problems. 

Two  types  of  adaptive  techniques  are  currently  popular:  (i)  moving  mesh 
methods,  where  a  grid  of  a  fixed  number  of  finite  difference  cells  or  finite 
elements  is  moved  so  as  to  follow  and  resolve  local  nonuniformities  in  the 
solution,  and  (ii)  local  refinement  methods,  where  uniform  fine  grids  are 
added  to  coarser  grids  in  regions  where  the  solution  is  not  adequately 
resolved.  A  representative  sample  of  both  types  of  methods  is  contained  in 
Babuska,  Chandra,  and  Flaherty  [1].  Roughly  speaking,  moving  mesh  methods  are 
superior  at  reducing  dispersive  errors  in  the  vicinity  of  wave  fronts  while 
local  refinement  methods  can,  in  principle,  add  enough  fine  grids  to  resolve 
any  fine  scale  structure  (cf.  Hedstrom  and  Rodrigue  [9]). 

We  discuss  a  moving  mesh  finite  element  procedure  for  finding  numerical 
solutions  of  m-dimensional  vector  systems  of  partial  differential  equations 
having  the  form 


(1.1)  Lu  :=  M(x,t)ut  +  f(x,t,u,ux)  -  [D(x, t,u)ux]x  *  0, 


0  <  x  <  1 ,  t>0. 


subject  to  the  initial  and  boundary  conditions 


t-'  < 


(1.2) 


u(x,0)  ■  u®(x)  ,  0  <  x  <  1  , 


( 1 . 2b,c)  either  u^x.t)  -  aj.(t)  or  £  Dj4U.s  (x,t)  -  a^Ct) 

j»l  x 

at  x  -  0  and  1  ,  t  >  0  ,  i  ■  1,  2,  ...,  m. 


We  are  primarily  concerned  with  parabolic  problems  where  D  and  M  are  positive 
definite;  however,  we  do  not  restrict  ourselves  to  this  case,  but  instead  we 
assume  that  conditions  are  specified  so  that  equations  (1.1,2)  have  an 
isolated  solution. 

We  discretize  (1.1,2)  in  space  using  a  finite  element -Galerkin 
procedure  with  piecewise  linear  approximations  on  a  moving  mesh.  We  simulta¬ 
neously  calculate  an  error  estimation  using  a  piecewise  quadratic  correction. 
This  error  estimate  is  used  to  move  the  mesh  so  that  it  approximately 
equidis tributes  the  local  spatial  component  of  the  discretization  error. 
Temporal  integration  is  performed  using  a  code  for  stiff  ordinary  differential 
equations  and  algebraic  systems  due  to  Petzold  [15].  We  halt  the  temporal 
Integration  at  specified  times  and  examine  our  error  estimate.  If  it  is 
larger  than  a  prescribed  tolerance,  the  step  is  rejected  and  the  integration 
is  re-done  using  a  finer  spatial  discretization.  On  the  other  hand,  if  the 
error  estimate  indicates  that  the  solution  is  being  calculated  too 
accuarately,  then  the  integration  is  continued  with  a  coarser  spatial  mesh. 

Our  procedure  differs  from  the  moving  finite  element  method  of  Miller  et 
al.  [8,12,13]  in  that  we  move  the  mesh  so  that  the  spatial  error  in  Hi  is 
equidistributed  and  they  move  their  mesh  so  as  to  minimize  the  residual  in 
1*2.  Our  refinement  procedure  also  differs  from  local  refinement  procedures 
of,  e.g.,  Berger  [2],  Bieterman  and  Babuska  [3,4],  and  Flaherty  and  Moore 
[6].  Our  mesh  equidistributes  the  local  error  and,  thus,  when  refinement  is 
necessary  it  is  performed  globally.  This  avoids  the  use  of  complicated  tree 
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data  structures  that  are  necessary  with  the  local  refinement  schemes  of,  e.g., 
Berger  [21  and  Flaherty  and  Moore  [6,71. 


In  Section  2  of  this  paper  we  discuss  our  discretization  procedure,  in 
Section  3  we  describe  our  algorithm  and  mesh  refinement  techniques,  in  Section 
4  we  apply  a  code  based  on  our  algorithm  to  several  linear  and  nonlinear 


examples,  and  in  Section  5  we  discuss  our  results  and  suggest  some  future 
considerations  and  refinements. 

2 .  Discrete  Formulation 

We  construct  a  weak  form  of  (1.1,2)  in  the  usual  manner.  Thus,  we  assume 

u  C  H^,  where  the  subscript  E  denotes  that  u  satisfies  any  essential 
E 

(Dirichlet)  boundary  conditions  in  (1.2).  We  select  a  test  function  vC  H*, 

0 

where  the  subscript  0  indicates  that  v  satisfies  homogeneous  versions  of  any 
essential  boundary  conditions.  We  multiply  (1.1)  by  v,  integrate  it  on 


0  <  x  <  1,  and  then  integrate  the  diffusive  terms  by  parts  to  obtain 


Accession  For 


(2.1a)  (v,ut)  +  ( v , f )  +  a(v,u)  -  0  ,  for  all  v  G  H1  ,  t  >  0  ,  KTTS  CRA&I 

n  r  t  a  a, 


where 


(2.1b) 


(2.1c) 


(v,u)  -  /  v(x,t)Tu(x,t)dx  , 


1  1 

a(v,u)  *  /  vT  D(x,t,u)uxdx  -  vTD(x,t,u)ux  | 
Ox  0 


r . iC  TAB 
C:  ii:i-i.;unced 

.Ti.ii  .Tlcntlor 


0\  •'  ••  i  S~ ■  it  i  on/ 

‘  •  ■  1  >1' lity  Codes 
Av».  11  aao./or 
3 1.>  <=•<.' ini 


To  construct  finite  element  solutions  of  (2.1)  we  select  finite 

dimensional  approximations  U  CS1*  and  V  C  to  u  and  v,  respectively. 

E  0 

Here,  SN  and  SN  are  finite  dimensional  subspaces  of  H*  and  , 

E  0  E  0 

respectively.  Thus,  we  seek  to  find  U  C  SN  satisfying 


\\Y 


4 

(2.2a)  (V,Ut)  +  (V,f)  +  a(V,U)  -  0  ,  for  all  V  CsJ  ,  t  >  0  , 

(2.2b)  (V,U)  -  (V,uO)  ,  t  «  0  . 

We  introduce  a  partition 

(2.3)  JI(t,N)  {0  -  xq( t )  <  xj(t)  <  ...  <  xN(t)  -  1} 

of  (0,1)  into  N  subintervals  (x^_j(t) ,x^(t)) ,  i  ■  1,  2,  ...»  N  ,  t  >  0  and 
select  U  and  V  to  be  piecewise  linear  polynomials  with  respect  to  this 
partition.  Thus,  for  example,  U  has  the  form 


N 

(2.4a)  U(x,t)  -  l  U1(t)4.i(x,n(t,N))  , 
i-0 

where 

f  *  -  x1-1(t) 

Xj[  ( t )  -  Xi-l(t) 


xi_1(t)  <  x  <  xi(t) 


(2.4b)  4^(x,n(t,N)) 


xi+l(t)  -  x 
xj.+i(t)  -  xt(t) 


x±(t)  <  x  <  *i+1(t) 


0  ,  otherwise 


We  note  that  the  derivative  Uc  in  (2.2a)  has  the  form 


(2.5) 


U 

t 


N  .  i+1 

I  {U.(t>*  (x,il(t,N))  +  l  V  (t) 
i-0  1  1  j-i-1  1 


<*<Ei 

dXj 


Xj(t)}  , 


where  (  )  :■  d(  )/dt.  Thus,  if  the  mesh  positions  and  velocities,  x^(t)  and 
Xj/t)  were  known,  we  would  have  a  set  of  ordinary  differencial  equations  for 
the  nodal  values  U<(t)  of  U(x,t). 


In  order  to  estimate  the  error  in  the  piecewise  linear  finite  element 
equations  (2.2)  to  (2.5)  we  write 


(2.6) 


u(x,t)  -  U(x,t)  +  e(x,t) 


and  substitue  this  into  (2.1a)  to  obtain 


(2.7)  (v,Ut+  et)  +  (v,f(.,t,U+e,Ux  +  ex))  +  a(v,0+e)  *■  0  , 

for  all  v  C  H1  ,  t  >  0  . 

0 

We  select  a  finite  dimensional  approximation  E  C  of  e  consisting  of 
piecewise  hierarchic  quadratic  functions,  i.e.. 


(2.8a) 


where 


E(x,t)  -  l  E1(t)ii(x,II(t,N))  , 
i-1 


(2.8b)  T^(x,II(t,N)) 


4[x  -  xi,i(t)}[xt(t)  -  x] 
{ x± ( t )  -  xi-i(t)Ji 


Xi-j(t)  <  x  <  xj.(t) 
otherwise 


We  also  select  a  finite  dimensional  approximation  V  G  of  v  consisting  of 
piecewise  quadratics.  Thus,  we  solve 


(2.9a)  (V,Ut  +  Et)  +  (V,f(.,t,U+E,Ux  +  Ex))  +  a(V,U+E)  -  0  , 


for  all  V  C  SN  ,  t  >  0  , 
0 


subject  to  the  initial  conditions 


(2.9b) 


(V,E)  -  (V,uO-U)  ,  for  all  V  C  S1*  ,  t  -  0  . 


Once  again,  if  the  mesh  positions  and  velocities  were  known,  we  would  have 
a  set  of  ordinary  differential  equations  for  E^(t).  Note  that  our  error 
estimate  E(x,t)  assumes  the  superconvergence  of  the  piecewise  linear  finite 
element  solution  U(x,t),  i.e.,  we  assume  that  the  piecewise  linear  solution 
is  converging  to  higher  order  at  the  nodal  points  x^(t),  i  -  0,  1,  ...,  N. 

The  error  estimate  that  is  calculated  by  solving  (2.9)  is  used  to 
control  the  motion  of  the  mesh.  Specifically,  we  determine  the  mesh  positions 
by  solving  the  ordinary  differential  system 

(2.10a)  x±(t)  -  Xi-i(t)  -  -X(|  |  Ei?ilh  -  E)  ,  i  -  1,  2,  ...  N, 

where  X  is  a  positive  constant  and 

_2  2  12  2 

(2.10b)  NE  -  I!  E  ||  :-  /  1(E)  +  (Ex)  ]dx  . 

Thus,  ||Ei¥i||i  is  the  local  error  in  H*  on  (x£_j(t),  x^(t))  and  E  is 
the  average  error  in  H1 .  If  | | Ei | | i  is  larger  than  E  then  the  right  hand 
hand  side  of  (2.10a)  is  negative  and  the  points  x^(t)  and  xi-l(t)  are 
moved  closer  together.  Similarly,  when  | | | | i  is  smaller  than  E  the  points 
x^(t)  and  X£_i(t)  are  moved  apart.  The  differential  system  (2.10)  was 
studied  by  Coyle  et  al.  [5]  and  shown  to  asymptotically  equidis tribute  the 
quantity  and  to  be  stable  relative  to  small  perturbations  in  the 

mesh  positions.  Large  values  of  the  parameter  X  give  shorter  relaxation 
times  of  x^(t),  i  *  0,  1,  ...,  N,  to  an  equidis tributing  mesh;  however, 
they  introduce  stiffness  into  the  system  (2.10).  This  strategy  is  similar  to 
one  suggested  by  Hyman  and  Naughton  (10]  and  we,  like  they,  solve  it  by 
eliminating  E  using  (2.10)  on  two  neighboring  intervals.  This  gives 


(2.11)  xi+1  -  2xi  +  xj.-!  -  -A(|  lEi+l’J'i+il  ll  ~  1 1  Ei  Vi  1  1 1 )  , 

i  -  1,  2,  ....  N-l  ,  t  >  0  . 

We  select  an  Initial  mesh  for  (2.11)  that  is  either  uniform  or  that 
equidistributes  E(x,0).  An  algorithm  for  calculating  an  equidistributing  mesh 
is  described  in,  e.g.,  Coyle  et  al.  [5]. 

At  present,  we  solve  the  three  sets  of  ordinary  differential  equations 
(2.2),  (2.9),  and  (2.11),  respectively,  for  the  piecewise  linear 
approximation,  the  quadratic  error  estimate,  and  the  mesh  positions 
simulataneously  using  the  backward  difference  code  DDASSL  developed  by  Petzold 
[15].  This  code  is  capable  of  solving  systems  of  implicit  differential  and 
algebraic  equations  having  the  general  form 

(2.12)  g(t,  y(t) ,  y(t))  -  0  ,  t  >  0  , 

and  subject  to  appropriate  initial  conditions.  It  can  solve  stiff  systems, 
with  variable  order  of  accuracy  and  temporal  step  size,  and  can  even  solve 
problems  when  the  matrix  M  in  (1.1)  is  singular.  Our  differential  system  is 
simpler  than  (2.12)  and  has  the  general  form 

(2.13)  A(t,y)y  +  h(t,y)  -  0  ,  t  >  0  . 

Furthermore,  the  matrix  A  and  the  Jacobian  of  h  are  sparse  and  banded  and 

DDASSL  is  capable  of  exploiting  this  structure.  At  present,  we  use  finite 

• 

difference  formulae  to  approximate  Jacobians  of  g  with  respect  to  y  and  y. 

3.  Mesh  Refinement  Algorithm 


set  of  output  times  toutput[k],  k  *  0,  1,  . 
tol,  and  a  value  for  X. 


K,  a  spatial  error  tolerance 


PROCEDURE  mferef (toutput,  K,  tol,  X); 

BEGIN 

select  initial  mesh; 
k  :=  1; 

WHILE  k  <  K  DO 
BEGIN 

solve  the  problem  on  toutput [k-1]  to  touput[k]  using 
the  moving  finite  element  procedure  described  in 
Section  2; 

compute  E; 


IF  E  >  tol 

THEN  add  a  mesh  point  to  each  subinterval 
ELSE 
BEGIN 


END 

END; 


IF  E  <  tol/10  THEN  delete  a  mesh  point  from 
each  subinterval; 


k  k  +  1 


END 


Figure  1.  Top  level  Algorithm  for  the  moving  finite  element 
procedure  with  mesh  refinement. 


We  halt  the  temporal  integration  at  every  output  point  and  examine  the 
average  error  E.  If  E  is  larger  than  tol,  we  add  a  finite  element  uniformly 
to  each  subinterval  and  re-do  the  integration.  If  E  <  tol/10  we  delete  every 
other  element  and  continue  the  integration.  In  the  case  when  tol/10  <  E  <  tol 
we  continue  the  integration  with  the  same  number  of  elements. 

Initial  conditions  are  needed  whenever  elements  have  been  added  or 
deleted.  When  refinement  is  necessary,  we  calculate  a  refined  solution  Ur(x,t) 
by  interpolating  U  +  E  at  time  t  using  (2.4a)  and  (2.8a),  l.e.. 


N  N 

(3.1a)  Ur(x,t)  -  l  UiCO^Cx.Itft.N))  +  l  Ei(t)  ^(x.lKt.N))  . 

i=0  i=l 

If  one  nodal  point  is  added  to  each  subinterval,  then  we  need  values  of 
Ur(xi  ,t),  i  *  0,  1,  N  and  Ur((xi  +  Xi+i)/2,t),  i  *  0,  1,  N-l. 

We  also  need  a  refined  error  estimate  Er(x,t)  and  we  calculate  this  as 

N  N 

(3.1b)  Er(x,t)  =  l  Et(t)  fi(x,  Il(t ,N) )  -  l  Ei(t)4>2i-l(x,n(t,2N))  . 

i=l  i=l 

Again,  if  one  nodal  point  Is  added  to  each  subinterval,  then  we  need  values  of 

Er((xi+Xi+i)/4,t)  and  Er(3(xi+xi+l)/4,t),  i  *  1,  2 . N. 

When  mesh  points  are  deleted,  we  need  values  of  E  on  the  coarser  mesh. 

If  every  other  mesh  point  is  deleted  then  we  use 

(3.2)  E(x2i_i,t)  -  U(x2i-i,t)  -  (U(x2i,t)  +  U(x2i-2tt) 1/2, 

i  -  1,  2 . N-l. 

Initial  values  for  U^(Q)  and  Ei(0)  were  obtained  by  interpolating  the 
exact  initial  function  u0(x).  The  initial  mesh  11(0,  N)  can  be  selected  so  as 
to  equidis tribute  E(x,0).  Finally,  we  select  a  constant  value  of  A  in  a 
relatively  ad  hoc  manner;  however,  we  are  studying  ways  of  automatically 
choosing  this  parameter. 

Procedures  for  evaluating  f(x,t,u,ux),  M(x,t),  D(x,t,u),  and  the  initial 
and  boundary  conditions  must  also  be  provided.  A  temporal  error  tolerance  and 
other  parameters  that  are  required  by  the  code  DDASSL  are  set  internally  in  a 
relatively  problem  independent  manner. 

4 .  Examples 


We  used  a  code  that  is  based  on  the  algorithm  of  Section  3  to  solve  the 
following  three  examples  that  Illustrate  the  performance  of  the  moving  finite 


element  method  and  the  utility  of  the  error  estimation  strategy.  In  order  to 
estimate  convergence  rates,  we  solved  some  of  the  examples  without  refinement 
and  on  stationary  meshes. 

Example  1.  We  consider  the  simple  linear  heat  conduction  problem 

(4.1)  ut  +  ux  -  uxx  *  f(x,t)  ,t>0  0<x<  1.0  , 

and  select  the  source  term  f(x,t),  the  initial  function  u0(x),  and  Dirichlet 
boundary  conditions  so  that  the  exact  solution  of  (4.1)  is 

(4.2)  u(x,t)  *  (1  -tanh(Ci(x-C2t-C3) 1 }/2  . 

The  solution  (4.2)  is  a  traveling  wave  and  its  steepness,  speed,  and 
phase  can  be  determined  by  selecting  C|,  C£,  and  C3,  respectively. 

We  solved  this  problem  on  stationary  and  moving  meshes  with  a  fixed 
number  of  elements.  In  each  case  we  chose  Cj  *  10,  C2  ■  1,  C3  *  -0.15,  and, 
for  the  moving  mesh,  A  *  10.  Our  results  for  the  exact  error  j|e(.,t)||},  and 
the  effectivity  Index 

(4.3)  0  -  |jE(.,t)||1  /||e(.,t)||1 

at  t  »  0.5  are  presented  as  functions  of  the  number  of  elements  N  in  Table  1. 
The  mesh  trajectories  for  N  **  20  are  shown  for  0  <  t  <  1.4  in  Figure  2.  An 
examination  of  Table  1  shows  that  )|e]J|  *  0(1/N)  for  both  fixed  and  moving 
meshes,  which  is  as  expected  [16].  The  error  of  the  moving  mesh  solution  is 
about  half  of  the  error  of  the  fixed  mesh  solution  when  N  <  40.  Furthermore, 
the  effectivity  index  0  ♦  1  as  N  increases  and  this  indicates  that  the 
error  estimate  E  converges  to  the  true  error  in  Hi  .  The  mesh  is  concentrated 
in  the  vicinity  of  the  wave  front  and  follows  the  wave  with  approximately  the 


correct  speed. 


We  next  solve  this  problem  on  a  moving  mesh  with  refinement.  We  select 
an  error  tolerance  of  0.1,  X  *  30,  and  show  the  solution  U(x,t)  and  the  mesh 
trajectories  in  Figures  3  and  4,  respectively.  Elements  are  added  as  the  pulse 
enters  the  region  0  <  x  <  1  for  small  t  and  then  they  are  deleted  as  the 
pulse  leaves  the  region  when  t  Is  about  1.25. 

Example  2.  We  consider  the  following  problem  for  Burgers'  equation: 
uc  +  uux  -  euxx  *0,  t>0,  0  <  x  <  1  . 

(4.4) 

u(x,0)  =  u°(x)  ,  u(0,t)  =  a(t)  ,  u(l,t)  *  b(t)  . 

As  a  first  example,  we  select  u^,  a,  and  b  so  that  the  exact  solution  of  (4.4) 
is  the  travleling  wave 

(4.5)  u(x,t)  *1-2  /e  tanhl (x-t)//e]  . 

As  in  Example  1,  we  solved  this  problem  on  stationary  and  moving  meshes 
with  a  fixed  number  of  elements.  We  chose  e  *  0.01  and,  for  the  moving  mesh, 

X  *  10.  Our  results  for  the  exact  error  J } e( . , t ) | | ^  ,  and  the  effectivity 
index  at  t  *  0.5  are  presented  as  functions  of  the  number  of  elements  N  in 
Table  2.  We  see  that  ||e||j  *  0(1/N)  for  both  fixed  and  moving  meshes  and 
that  the  error  estimate  converges  to  the  true  error.  Here,  as  in  Example  1,  we 
see  that  the  effectivity  index  is  converging  at  a  faster  rate  for  the  uniform 
mesh  than  for  the  moving  mesh. 

As  a  second  example  for  (4.4),  we  select 

(4.6)  u0  (x)  -  10x(x  -  1 ) (x  -  3/4)  ,  a(t)  -  b(t)  -  0  . 

For  small  times  and  e,  the  exact  solution  is  a  pulse  that  moves  in  the 
positive  x  direction  while  steepening.  At  about  t  ■  1,  a  shock  layer  forms 
near  x  *  1  and  after  a  time  of  0(l/e),  the  solution  dissipates  to  zero. 


STATIONARY  MESH 

MOVING 

MESH 

N 

mimm 

e 

mam 

0 

10 

0.4275 

0.9634 

0.2810 

0.9232 

20 

0.2330 

0.9941 

0.1604 

0.9874 

40 

0.1175 

0.9985 

0.0892 

0.9978 

TABLE  1.  Exact  error  ||e||i  and  effectively  index  0  for  Example  1 
Calculations  well  performed  on  stationary  and  moving  meshes  having 
N  elements. 


STATIONARY  MESH 

MOVING 

MESH 

N 

mm 

6 

1 lel ll 

0 

10 

0.2036 

0.8027 

0.1499 

0.7321 

20 

0.09372 

0.9724 

0.0847 

0.9144 

40 

0.04707 

0.9950 

0.0434 

0.9854 

TABLE  2.  Exact  error  | | e | [ ^  and  effectively  index  0  for  Example  2 
with  the  exact  solution  given  by  eq.  (4.5).  Calculations  were 
performed  on  stationary  and  moving  meshes  with  N  elements. 


We  preseat  the  solution  U(x,t)  and  the  mesh  trajectories  for 
computations  performed  with  20  moving  elements ,  e  *  0.01 ,  and  X  *  10  in 
Figures  5  and  6,  respectively.  The  mesh  is  concentrated  in  regions  of  high 


curvature  and  follows  the  moving  wave  front. 

Example  3.  We  consider  the  reaction-diffusion  model  in  one  dimension 
(cf.  Kapila  [11]) 

ut  «  uxx  -  Due-^/^  , 

(4.7a,b)  t>0,0<x<l  . 

LTt  "  Txx  +  aDu  e“6/T  , 

e6 

(4.7c)  D  -  R 


(4.7d,e)  u(x,0)  -  T(x,0)  «  1  , 

(4.7f,g)  Uj  (0,t)  -  Tx(0,t)  »  0  , 

(4.7h,i)  u(l,t)  -  T( 1 ,t)  -  1  . 

This  model  discribes  a  single  one  step  reaction  (A  +  B)  of  a  gas  in  a 
region  0  <  x  <  1.  The  quantity  u  is  the  mass  fraction  of  the  reacting  gas,  T 
is  the  gas  temperature,  L  is  the  Lewis  cumber,  a  is  the  heat  release, 

4  is  the  activation  energy,  D  is  the  Damkohler  cumber,  and  R  >  0.88  is  the 
reaction  rate. 

We  first  consider  th%  case  when  L  ■  1.  Then  T  +  cu  *  1  +  a  and 
(4.7)  reduces  to  the  simpler  problem 


(4.8a) 


Tt  *  Txx  +  D(l+cr-T)e-S/T  ,t>0,0<x<l. 


(4.8b,c)  Tx(0,t)  -  0  ,  T( 1 ,t)  -  1  ,  t  >  0  , 

■  l  .  n  t.  v  c  i  . 


For  small  times  the  temperature  gradually  increases  from  unity  with  a  "hot 
spot"  forming  at  x  *  0.  At  a  finite  time,  ignition  occurs  and  the  temperature 
at  x  *  0  jumps  rapidly  from  near  unity  to  1  +  a.  A  sharp  flame  front  then 
forms  and  propagates  towards  x  *  1  with  speed  proportional  to  e°^/2(l+a).  In 
real  problems,  a  is  about  unity  and  £  is  large;  thus,  the  flame  front  moves 
exponentially  fast  after  ignition.  The  problem  reaches  a  steady  state  once  the 
flame  propagates  to  x  *  1. 

We  solve  (4. 7c, 8)  for  a  “  1,  6  *  20,  and  R  *  5  using  a  moving 
mesh  with  20  elements.  Our  results  for  the  mesh  trajectories  when  X  *  10  and 
200  are  shown  in  Figures  7  and  8,  respectively.  The  computed  temperatures  T 
vs.  x  are  shown  in  Figure  9  for  several  times  and  X  *  10  and  200.  We  also 
computed  solutions  on  a  moving  mesh  with  refinement  using  tol  “0.1  and 

X  “  300.  The  mesh  trajectories  and  solution  in  this  case  are  shown  in  Figures 

10  and  11,  respectively. 

We  see  that  the  ignition  time  and  the  computed  solutions  while  the  flame 
front  is  propagating  from  x  «  0  to  x  *  1  are  sensitive  to  the  choice  of  X  (cf. 
Figures  9  and  11).  This  is  due  to  the  extremely  fast  rate  of  ignition  and 
velocity  of  the  front.  All  three  solutions  are  in  essential  agreement  before 
ignition  and  upon  approaching  steady  state.  This  is  a  very  difficult  problem 
and  yet  our  methods  are  capable  of  finding  a  solution  with  relative  ease.  This 
is  because  the  mesh  is  moving  with  approximately  the  speed  of  the  front  and 
the  solution  along  the  mesh  trajectories  is  changing  slowly. 

As  a  final  problem,  we  solve  the  coupled  system  (4.7)  with  L  «  0.9.  As 
in  the  previous  example,  we  take  a  «  1,  6  “  20,  R“5  and  solve  the  problem  on 

a  moving  mesh  with  20  elements.  Our  results  for  the  mesh  trajectories  when 

X  “  10  and  150  are  shown  in  Figures  12  and  13,  respectively.  The  computed 


temperature  and  mass  fraction  are  shown  in  Figure  14  for  X  -  10  and  in  Figure  15 

for  X  =  150.  Again  we  computed  a  solution  with  refinement  using  tol  ■  0.1  and  X 

»  500.  The  mesh  trajectories  and  solution  are  shown  in  Figures  16  and  17, 
respectively.  Once  again,  the  ignition  time  and  computed  solutions  have 
different  phases  while  the  flame  front  is  moving.  We  show  the  estimated  error 
| |e| ji  as  a  function  of  time  for  solutions  calculated  with  X  *  10,  N  -  20,  X  ■ 

150,  N  «  20,  and  X  *  500  with  refinement  in  Figure  18.  For  N  ■  20,  the  estimated 

error  is  about  50  percent  smaller  when  X  =  150  than  for  X  «*  10.  The  estimated 
error  is  less  than  the  prescribed  tolerance  of  0.1  when  refinement  is  used. 

5 .  Discussion  of  Results 

We  have  presented  a  moving  finite  element  method  where  the  mesh  is  moved 
so  as  to  equidis tribute  the  spatial  component  of  the  discretization  error  in 
H*.  We  also  discuss  a  computationally  simple  and  effective  procedure  for 
estimating  the  discretization  error  using  quadratic  hierarchic  finite 
elements.  The  error  estimate  gives  users  some  confidence  in  the  computed 
solution  and  we  also  use  it  to  develop  a  simple  procedure  for  combining  mesh 
movement  and  refinement.  Our  approach  makes  use  of  existing  software  for 
solving  ordinary  differential  equations .  We  have  used  the  differential 
algebraic  system  solver  DDASSL  [15J;  however,  there  are  several  other 
possibilities .  We  have  applied  a  code  based  on  our  method  to  several  examples 
and  have  shown  that  it  can  effectively  solve  some  very  difficult  partial 
differential  systems  with  little  or  no  user  intervention. 

While  our  results  are  very  encouraging,  there  are  several  aspects  of  our 
work  that  are  still  incomplete.  For  some  problems  (e.g.,  Example  3)  the 
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Figure  9.  Temperature  T(x 
for  Example  3, 

( diamonds ). 


X 

Figure  14.  Temperature  (triangles)  and  mass  fractioo  (octagons)  vs.  x  at 
t  -  0.0,  0.227,  0.229,  0.237,  0.245,  and  0.260  for  Example  3, 
Eq.  (4.7),  with  L  ■  0.9,  X  «  10  and  N  ■  20. 


X 

Figure  15.  Temperature  (triangles)  and  mass  fraction  (octagons)  vs.  x  at 
t  ■  0.0,  0.227,  0.229,  0.237,  0.245,  and  0.260  for  Example  3, 
Eq.  (4.7),  with  L  ■  0.9,  X  *  150  and  N  *  20. 


